Trapped fermions with density imbalance in the BEC limit 
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We analyze the effects of imbalancing the populations of two-component trapped fermions, in 
the BEC limit of the attractive interaction between different fermions. Starting from the gap 
equation with two fermionic chemical potentials, we derive a set of coupled equations that describe 
composite bosons and excess fermions. We include in these equations the processes leading to the 
correct dimer-dimer and dimer-fermion scattering lengths. The coupled equations are then solved 
in the Thomas-Fermi approximation to obtain the density profiles for composite bosons and excess 
fermions, which are relevant to the recent experiments with trapped fermionic atoms. 

PACS numbers: 03.75.Hh,03.75.Ss 
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Lately there has been mounting interest in trapped 
atomic fcrmion superfluids, both experimentally and the- 
oretically These systems allows one to explore the BCS- 
BEC crossover with the use of Fano-Feshbach resonances. 
Interesting many-body issues for fermionic superfluids 
have been addressed in this context, aiming also at mak- 
ing connections with high-temperature superconductors 
and nuclear and quark matters. Until recently, only the 
case of two fermionic species with equal populations was 
experimentally explored. 

Two new experimental studies of fermionic superfluid- 
ity with trapped atoms having imbalanced populations 
yj, |2| have raised novel interest in these systems. Den- 
sity profiles of the two fermionic species 1, 2] as well as 
vortices Q have been detected. A quantum phase tran- 
sition to the normal state on the weak-coupling side of 
the BCS-BEC crossover pj as well as a phase separa- 
tion in the crossover region P, Q , both driven by density 
imbalance, have been identified. 

The many-body problem becomes reacher when the 
fermionic populations are imbalanced, since this addi- 
tional degree of freedom gives rise to new phenomena. 
The effects of density imbalance on fermionic superfluids 
have been studied theoretically mostly within mean field 
and on the weak-coupling (BCS) side of the crossover, 
both for the homogeneous [3| and trapped case |4| . Only 
recently these calculations have been extended to cover 
the BCS-BEC crossover [j| , but without considering the 
effects of the trap which are essential to account for the 
experimental results with density imbalance |lj, |2j ■ 

The effects of the trap are especially important to pro- 
duce a phase separation between the density profiles of 
the two fermionic species, which may otherwise remain at 
the verge of an instability for a homogeneous system. In 
this respect, consideration of the strong-coupling (BEC) 
side of the crossover appears relevant, since phase separa- 
tion is more robust in this limit. [On the weak-coupling 
side, on the other hand, density imbalance acts to quickly 
destroy the superfluid phase.] In addition, on the BEC 
side theoretical studies can rely on more accurate treat- 
ments beyond mean field, by exploiting the dilutcness 



condition of the system. From previous experience on the 
BCS-BEC crossover, one expects the results obtained on 
the BEC side to be qualitatively similar to those occur- 
ring near unitarity. 

Accordingly, in this paper we analyze the BEC side of 
the crossover where trapped fermions with imbalanced 
populations get rearranged into a system of compos- 
ite bosons (dimers) plus a number of excess (unpaired) 
fermions. We consider specifically the limit of low tem- 
perature, whereby all dimers are condensed, and derive a 
set of coupled equations describing the dimers and the ex- 
cess fermions, which are interacting through dimer-dimer 
and dimer-fermion scattering processes. To the extent 
that the system is dilute, we can consider these processes 
exhaustively. Due to these interactions, the mutual ef- 
fects between dimers and excess fermions need to be dealt 
with self-consistently in the equations. This results in a 
nontrivial evolution of their density profiles as the degree 
of density imbalance increases. We determine numeri- 
cally these density profiles within a local approximation 
and find that dimers and excess fermions tend to reside 
in different spatial regions of the trap. We also show how 
from these profiles one can, in reverse, extract the ra- 
tio Obf /«f of the dimer-fermion to the fermion-fermion 
scattering lengths. 

We first consider the derivation of the above coupled 
equations within mean field. When the populations of 
the two fermionic species (labeled by | and J.) are equal, 
it has already been shown ||| that the Bogoliubov-de 
Gennes (BdG) equations for trapped fermions can be 
mapped in the BEC limit onto the Gross-Pitaevskii equa- 
tion for composite bosons, the only remnant of the origi- 
nal fermionic nature residing in the dimer-dimer scatter- 
ing length entering that equation. In the case of interest 
here of imbalanced populations, this mapping has to be 
reconsidered since composite bosons and excess fermions 
coexist and mutually interact. In the following, we shall 
recall schematically the derivation of Ref. 6], and indi- 
cate where and how it needs to be modified to account 
for population imbalance. 

The presence of two fermionic species with popula- 
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tions TVj and requires us to introduce two chemi- 
cal potentials ^ and (For definiteness, we assume 
> N±.) In analogy to Ref.[||, we consider the non- 
interacting Green's functions which satisfy the equation 
[iuj s - H a (r)] Go(t,t';u b \h <t ) = 6(r-r'), being subject 
to the trap potential V(r) entering the single-particle 
Hamiltonian TL a {v) = — V 2 /(2m) + V(r) — \i a . Here, 
a = (T) I); m 1S the fermion mass, and uj s = (2s + l)7r//3 
(s integer) is a fermionic Matsubara frequency (/? be- 
ing the inverse temperature). With these noninteracting 
Green's functions we form the matrix: 



Go(r,r';uj s ) 



g (r,r';u s \fj,f) 

-eo(r',r;-w s |//j.) 



The corresponding interacting single-particle Green's 
functions within mean field are obtained by solving the 
integral equation: 

g(r,r';u s ) = £? (r, r'; w s ) 

dr"go(r,r";a; s )B(r")e(r",r> s ) (1) 



where 



B(r) = 



A(r) 
A*(r) 



contains the local gap A(r) . Solution of the integral equa- 
tion is equivalent to the solution of the BdG equations 
in the trap in the presence of density imbalance. 

The elements (1, 1) and (2, 2) of the matrix JJ) deter- 
mine the local densities, while the element (2,1) deter- 
mines the local gap. To obtain the coupled equations 
of interest we expand them up to order A 2 and A 3 , re- 
spectively, in analogy to what was done in Ref. 6] in the 
absence of density imbalance. The gap equation reads: 







- i& (r)V 2 A*(r)+ (- + a (r)) A*(r) 
+ c (r)|A(r)| 2 A*(r) 



(2) 



where r\ is a positive infinitesimal and vq (< 0) is the 
strength of the attractive contact potential between 
fermions of opposite spins. As in Ref.[6j, the coefficients 
ao, bo, and cq of Eq.JSJ are expressed in terms of the 
noninteracting Green's functions Go(fJ>a), here with chem- 
ical potentials fi a appropriate to imbalanced populations. 
While bo(r) = maF/(167r) retains the same value as for 
equal populations, ag and cq are affected by density im- 
balance in a relevant way. With the local ansatz for 
Go(^a) introduced in Ref.[|| in terms of the noninter- 
acting Green's functions of the associated homogeneous 
problem, we obtain at low temperature: 

1 1 f dk 

v v J (27r)' j 



2m 



V(r) 



m 2 ap 
8tt 



[Mi 



£ - MT - m + 2V(r) 
2V(r)] - —6n f (r) 



(3) 



where /(e) = (exp(/3e) + 1) _1 is the Fermi function. Wc 
have made use of the fermion bound-state equation with 
binding energy e = (ma|) _1 , introduced the bosonic 
chemical potential /ie = /i| + /i| + eo, and anticipated 
that /i| ~ — eo while /.if ~ ep(dn) where ep(i5n) is the 
Fermi level of free excess fermions with density: 



2m 



+ V(r) - M 



(4) 



Equation Q can be cast in the meaningful form 

MB _ 2 y (r) _ 3 ^I 5n/(r)i 

'(5) 

where we have identified the dimer-fermion scattering 
length cl-qf which takes the Born value 8ap/3 at the 
mean-field level. Wc obtain in addition: 



1 . . m 2 aF 
— +a Q r = — — 

V 87T 



co(r) 



_d d_ 



oo (r) = - 



77i dp 2 4 88nf(r) 



16tt 



+ m a-p 



<9/j t 



(6) 

The gap equation J2J eventually reduces to the follow- 
ing equation for the condensate wave function $(r) = 
y/m 2 a F /(8Tr) A(r): 

V 2 



2m B 



$(r) 



2T/(r) 
47raB 



37raBF 



Sn(r) 



$(r) 



777B 



$(r)| 2 $(r) =MB$(r) (7) 



where me = 2m is the dimer mass and the dimer-dimer 
scattering length takes the Born value 2of at the 
mean- field level. Note how the last contribution in Eq. JfjJ) 
combines with the term proportional to Sn f in Eq. JSJ , to 
yield in Eq.10) the density of excess fermions: 



5n(r) = 5rif(r) |<&(r)' 2 



(8) 



dk /k 2 

(2tt) 3 ; ^2771 



37ra B F 



|$(r)| 2 -/i T 



By this mechanism, Eq.Ql contains the effect of the 
fermionic distribution Sn(r) on the bosonic one, an effect 
which is reciprocated in Eq.© by the bosonic distribu- 
tion <£>(r)| 2 on the fermionic one. 

The final equation results from considering the total 
density n(r) — n^(r) + nj.(r) where: 



n ( r ) = pY, 

= 8n(r) - 



^'Gii(r,r;a; s 
2|$(r)| 2 . 



-^"G 22 (r,r;^ s )] 
(9) 
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FIG. 1: Tadpole diagrams for composite bosons representing 
the gap equation in the BEC limit: Mean-field contributions 
(a); Contributions including the dimer-fermion (b) and dimer- 
dimer (c) scattering processed beyond the Born approxima- 
tion. Full line: fermionic propagator of given spin; broken 
line: fermionic interaction vo- 



Equations J7J and JjJJ , together with the definition JHJ , 
are the desired closed set of equations which determine 
<&(r) and Sn(r) in the BEC limit, for given total number 
of fermions + N± = Jdrn(r) and density imbalance 
-/Vf — Ni = JdrSn{r). 

The validity of these equations can be extended by 
improving on the values of the scattering lengths obf 
and ae- At the mean-field level so far considered, these 
values correspond to the Born approximation, yielding 
&bf = 8ap/3 and ae = 2ap. The exact value 1.18of of 
czbf has been known for some time while the exact 
value 0.6ap of ob was determined only recently 8]. These 
values have also been obtained by diagrammatic methods 
in the limit of vanishing density, for obf hi Ref . |IJ and for 
<2b in Rcf. 10]. To improve on the derivation of Eqs.J7|) 
and @, so as to include the exact values of obf and ob, 
we thus need to identify additional fermionic diagrams 
beyond mean field, which contain the correct diagram- 
matic sequences for cibf and as as sub-units. This can 
be achieved as follows. 

In the BEC limit, the gap equation can be inter- 



preted as the condition of vanishing "tadpole" insertions 
for composite bosons, with a natural extension of what 
is done for point-like bosons [ll|. Within mean field, 
the diagrams representing this condition are depicted in 
Fig. 1(a). Here, a composite-boson propagator with zero 
four-momentum can be inserted from the left, while the 
gap A corresponds to a condensate line. We have already 
shown that these diagrams contain the dimer-fermion 
and dimer-dimer scattering processes within the Born 
approximation. Additional tadpole diagrams, whose ef- 
fect is to modify the values of obf and as in Eq. , are 
depicted in Figs. 1(b) and 1(c), in the order. These di- 
agrams exclude the Born contributions already included 
in the mean-field diagrams of Fig. 1(a). A comment is in 
order on how to factor out from the diagrams of Fig. 1(b) 
the product (Zbf times Sn that enters Eq.J7J. From the 
structure of these diagrams one concludes that the in- 
tegration over the wave vector P is bounded within the 
Fermi sphere with radius ^/2m^,-\, while the remaining 
integrations over the wave vectors q, q', • • •, extend out- 
side this Fermi sphere. Accordingly, one may neglect the 
P-dependence everywhere in the diagrams of Fig. 1(b), 
except in the fermion propagator labeled by P and cor- 
responding to a spin- 1 fermion. The density of excess 
fermions results in this way from the P-integration [l2j . 
while the remaining parts of the diagrams yield the exact 
dimer-fermion scattering matrix T 3 (0), which excludes 
the Born contribution resulting from mean field. 

A similar analysis can be carried out for the density 
equation J^J , where only the value of the dimer-fermion 
scattering length csbf needs to be corrected beyond mean 
field [cf. Eq.©]. In this case the relevant diagrams are 
similar to those of Fig. 1(b), but with an additional con- 
densate (A) line in the place of the incoming composite- 
boson propagator with zero four-momentum and with 
a doubling of the spin-| fermion propagator with four- 
momentum P. This doubling provides the derivative 
with respect to /if in Eq.JSJ, while the remaining parts 
of the diagrams yield again T 3 (0). 

We thus solve the coupled equations J7J and © with 
the exact values of obf and ob- To make contact with 
the recent experimental findings 0, 0] , we determine the 
density profiles Sn(r) and n (r) = $(r)| 2 by solving 
these equations in the Thomas-Fermi approximation for 
a spherical trap, as functions of the asymmetry parameter 
a = (iV T - 2V;)/(iV T + Ni) (with < a < 1). 

In Fig. 2 we report the radial density profiles Sn(r) and 
no(r) vs the distance r = |r| from the center of the trap, 
for three characteristic values of a and for the coupling 
(^fOf) 1 = 3 on the BEC side of the crossover. The spa- 
tial separation between the condensed composite bosons 
and the excess fermions is evident from these plots. For 
the coupling here considered, this phase separation oc- 
curs for all values of a, because Sn(r) tends to set outside 
the region occupied by n (r). Note that, for each value 
of a, the maximum of Sn(r) occurs where no(r) vanishes. 
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FIG. 2: Radial density profiles Sn(r) (full lines) and no(r) 
(broken lines) vs r = |r| for (fepap) -1 = 3 and a = 0.2 (upper 
panel), a — 0.5 (middle panel), a — 0.8 (lower panel) [r is in 
units of Rtf and densities are in units of (iVf + JVj, ) /-Rtf] • 
The insets show the results for (fepap) -1 = 1 at the same 
values of a. Here, Rtf and fep are the TF radius and Fermi 
wave vector for equal populations. 
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FIG. 3: Column density profiles Sn(p) (full lines) and no(p) 
(broken lines) vs p for (fepap) -1 = 3 and a — 0.2 (upper 
panel), a = 0.5 (middle panel), a = 0.8 (lower panel) [p is in 
units of Rtf and densities are in units of (iVf + Ni)/ R TF ]. 
The insets show the results for (fepap) -1 = 1 at the same 
values of a. Here, Rtf and feF are the TF radius and Fermi 
wave vector for equal populations. 



Note also the progressive size shrinking of no(r) for in- 
creasing a, with a simultaneous penetration of 6n(r) to- 
ward the center of the trap. The insets show the results 
for the smaller coupling (fepOF) -1 = 1- The phase separa- 
tion between the two species appears sharper in this case, 
corresponding to enhanced effects of the dimer-fermion 
repulsion. For the column density profiles reported in 
Fig. 3 (obtained from the above radial profiles by integrat- 
ing over, say, the z axis) the phase separation appears less 
pronounced even for (/cf«f) -1 = 1> since the correspond- 
ing distribution J dz 5n{p, z) leaks toward p = where it 
acquires a finite value. In particular, for the coupling 
(kvap)^ 1 = 1 near the boundary of the BEC region our 
column density profiles are similar to those reported in 
Fig. 2 of Ref.|2]. To enhance the visibility of phase sepa- 
ration, it would therefore be preferable to consider radial 
rather than column density profiles. 

Our calculation does not reveal evidence for a criti- 
cal value a c below which the phase separation does not 
occur. This somewhat contrasts with the experimental 
claim made in Ref.0, where at unitarity a c was esti- 



mated to be about 0.10. This difference may either be 
due to the limited experimental resolution when reveal- 
ing small excess densities, or possibly to the vanishing of 
a c when the BEC region (considered in our calculation) 
is approached. 

In Fig. 4 the critical radius R c of the condensate den- 
sity is plotted vs a for different couplings, showing specif- 
ically how the condensate density shrinks for varying a 
and coupling. This quantity should have direct experi- 
mental access, because it identifies also the position of the 
maximum of the density of excess fermions. This maxi- 
mum value can be compared with the value of the excess 
density at the center of the trap, to determine how their 
ratio v depends on the dimcr-fcrmion scattering length 
czbf (assuming that the value of ob is independently de- 
termined). Column density profiles are found to have 
a more marked dependence on obf than radial density 
profiles. Accordingly, in the inset of Fig. 4 we report the 
corresponding ratio v vs obf for a = 0.5 and different 
couplings. The rather marked dependence on obf shown 
by this quantity should make it possible to extract from 
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FIG. 4: Critical radius R c of the condensate density vs the 
asymmetry parameter a (normalized to the value i? c (0) at 
a = 0) for the coupling values (<:f1f) _1 = (1,2,3) from bot- 
tom to top. The inset shows the ratio v between the column 
density of excess fermions at the center of the trap and its 
maximum value at R c vs bbf/sf, for a — 0.5 and the cou- 
pling values (fcFaF) -1 = (2, 3, 4) from bottom to top. 



the experimental data the expected value 1.18 of obf/of, 
using our plots for calibration. 

We mention, finally, that, to compare with Fig. 2D of 
Ref.Q, we have also calculated the ratio v for a = 0.57 
and various couplings down to (/cf^f) -1 = 1, and extrap- 
olated eventually the results toward the unitarity limit, 
obtaining for this ratio about 0.5 in agreement with the 
experimental value. 
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